Dynamical Phase Transitions In Driven Integrate- And-Fire Neurons 
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We explore the dynamics of an integrate-and-fire neuron with an oscillatory stimulus. The frus- 
tration due to the competition between the neuron's natural firing period and that of the oscillatory 
rhythm, leads to a rich structure of asymptotic phase locking patterns and ordering dynamics. The 
phase transitions between these states can be classified as either tangent or discontinuous bifurca- 
tions, each with its own characteristic scaling laws. The discontinuous bifurcations exhibit a new 
kind of phase transition that may be viewed as intermediate between continuous and first order, 
while tangent bifurcations behave like continuous transitions with a diverging coherence scale. 



Neurons in awake, behaving mammals receive complicated 
dendritic input currents and respond with highly irregular 
trains of action potentials. Unraveling the meaning of each 
neuron's train of spikes is a formidable challenge. A simple 
starting point is to characterize a neuron's activity in terms 
of its firing rate (which varies in time in behaving organ- 
isms) , or in terms of temporal correlations between its spike 
times and that of other neurons, either individually or collec- 
tively in a local rhythm. There is a long history of detecting 
rhythmic neural activity at various scales, from electroen- 
cephalography (EEG), to local field potentials, to oscillatory 
membrane currents stimulating individual pyramidal and in- 
terneurons in voltage clamp recordings. 

In this letter we present a biophysical approach to explore 
consequences of rhythmic inputs on the rate and timing of 
a model neuron's spikes. Our framework is directly relevant 
for describing the behavior of a neuron in a slice prepara- 
tion with blocked dendritic inputs and controlled injection 
of a simple stimulus current in a whole-cell patch-clamp set- 
ting. Our analysis may also shed light on the time scales and 
(transient) dynamical patterns in the spike trains that may 
develop for the more complicated stimuli neurons receive in 
vivo. Advocating a statistical-mechanical perspective, where 
minimal models yield insight into universal behavior of more 
realistic models, we consider the one-dimensional integrate- 
and-fire (IF) modelpr. An IF model neuron receiving a con- 
stant current stimulus has a constant firing rate. We consider 
the consequences of an additional small oscillatory input; this 
second, competing time scale introduces frustration that re- 
sults in pattern formation. We explore connections between 
these stable patterns and the description of critical phenom- 
ena associated with continuous phase transitions. 

The IF model describes the response of a cell's mem- 
brane potential v{t) to an infiux of current /. The volt- 
age evolves according to the differential equation T{dv/dt) = 
— (v — Veq) + RI , together with the condition that when v{t) 
reaches a threshold vth, an instantaneous action potential is 
generated and v{t) is reset to an equilibrium resting poten- 
tial Veq- We use units where R is the membrane resistance 
and the time constant r = RC is proportional to the mem- 
brane capacitance C . The parameters R, I, r, Veq and Vth 



are all time-independent. We set Aw — vth 



Jeq 



and assume 



throughout that A?; > 0. If we start with an initial condition 



J, then the potential v{t) will always reach thresh- 



stant rate, with interspike interval Tnat — —t log(l— Aw/i?/). 

Next we consider the effect of adding a secondary periodic 
stimulus current E cos tut, which introduces a competing time 
scale Tdrv = ^tt/uj. The governing equation is now 



dv 

T— = — (w — Veq) + RI + EcOSUt. 



(1) 



We view the rhythm E cos ujt as a perturbation of the con- 
stant current term, and in this spirit limit our attention to 
the case < < RI. The introduction of a second, compet- 
ing time scale leads to a loss of the simple periodic behavior 
of the original model. As we shall see, the model neuron typ- 
ically no longer has a constant interspike interval, and can 
exhibit both periodic and aperiodic firing patterns. 



We again start with v{to) 



(changing allows us 



in effect to adjust the relative phase of the cosine drive at 
the initial condition). The next spike time ti is the first 
solution to v(t) = Vffi with t > to; a solution exists provided 
RI + Ej ^{ujtY + 1 > Aw. The dependence of ii on Iq 
determines a return map F\ iterating the map F generates a 
spike train < ii < ^2 • • • via t„ = F{tn_i) — i^"(io). 
The map F is continuous when 



RI>E + Av 



(2) 



(this implies v > for v < Vth) and is discontinuous other- 
wise. The source of this discontinuity is illustrated in Fig.[T| 
which shows solutions to ([l]) with nearly identical initial con- 
ditions in the two different cases where condition ^ is sat- 
isfied (green) and not satisfied (red). When ^ fails ([T]) has 
a solution with ii — Q ai threshold, which causes a jump 
discontinuity in the return map F . 

We are interested in the evolution of the interspike inter- 
vals (ISIs) tn — tn-i for E > 0, especially whether these 




old provided RI > Au. In this case the neuron fires at a con- 
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FIG. 1: Solutions to ([T]) when ([2| holds (green) and fails (red). 
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FIG. 2: Average period in units of the drive period vs. the pa- 
rameter RI for drive amplitude E = O.lAn, r — 20ms and 
Tdrv = 35ms. The dotted curve is T„at/Tdrv, which is propor- 
tional to the period in the undriven case. 

intervals become periodic and if so, how rapidly such a pat- 
tern is established. The asymptotic dynamics are determined 
by the average interspike interval 



^=lim^ 

n — *oo 77, n— >oo Jl 



(3) 



i.e., the inverse of the neuron's firing rate. 

The map F satisfies the periodicity relation F{t + Tdrv) = 
F{t) + Tdrv reflecting the periodicity of the drive term. Con- 
sequently, from the theory of circle maps [2 [3l |4] , the limit 
defining Tave exists, is independent of the initial condition 
to and depends continuously on the parameters {RI, r, lu 
and E). Furthermore, the dimensionless ratio Tave/Tdrv is a 
rational number r = p/q ]& 



drv 



(4) 



for some t*; in other words t* is a fixed point of the map 
F'^it) — p Tdrv So the spike train beginning with to = t* 
satisfies tn+q = tn+ pTdrv and consequently the sequence of 
phases of i„ relative to Tdrv repeats every q firings. 

In Fig. |2]wc plot Tave/Tdrv as a function of the parame- 
ter RI, while keeping fixed Tdrv — 35 nis, t — 2Q ms and 
E = 0.1 Av. We divide the graph into two regions accord- 
ing to whether the return map is continuous (green) or dis- 
continuous (red). For a given number r, let RI^ and 
denote the minimum and maximum values of RI for which 
Tave/Tdrv = f- Rlr < ^-^r^ iff r is a rational number; in 
other words the plateaux in Fig. [2] correspond to rational 
multiples of the drive period. When r is an integer and the 
return map is continuous (green part), RI^ can be deter- 
mined algebraically; in this special case 



RI^ — 



^ p — 27rr/ur 



± 



E 



(5) 



and T„. 



Tnat at the midpoint of the plateau. The width 



of such a plateau is then proportional to E, the magnitude 
of the oscillatory drive. In all other cases, RI^ needs to be 
determined numerically. 
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FIG. 3: (a) Return maps associated with the tangent bifurcation 
at left edge of the r = 1 plateau; (b) Return maps associated with 
the discontinuous bifurcation at right edge of the r = 2 plateau. 

The asymptotic structure of the average firing rate as in 
Fig. |2]has been known for some time[3J IH |S]. In this paper 
our focus is on the approach to the asymptotic behavior, and 
the resulting connection to dynamical phase transitions. The 
circle map theorem [2] guarantees that the asymptotic value 
of Tave for a given RI is independent of the initial condition 
to- This stability allows us to view Fig.[2]as a phase diagram, 
with each plateau a state corresponding to some rational 
number r which we express as a fraction r = p/q (in lowest 
terms). These states are analogous to phases of matter and 
the boundaries of the plateaux to phase transitions. 

Within each p/q entrainment plateau, the spike train con- 
verges to a periodic pattern of ISI's that repeats every q 
spikes, corresponding to a stable fixed point of F'^(t) —pTdrv 
For perfect p/q entrainment tn+q ~ tn ~ pTdrv — 0; hence 



AP'« = t 



n+q 



drv 



(6) 

measures the deviation from p/q entrainment. The conver- 
gence (within a plateau) is geometric, so A^'"^ ~ x~^. 

Henceforth we fix r, uj and E and consider the parametric 
dependence on RI. The fixed points of F'^{t) ~ pTdrv vary 
with RI and ultimately vanish through some kind of bifurca- 
tion at the edges of the p/q entrainment plateau. Analogous 
to the theory of phase transitions, the equation 



(7) 



defines a coherence time ^r{RI) that characterizes how 
rapidly the phase-locked solution is approached. Of particu- 
lar interest is how ^^(i?/) scales with respect to the tuning 
parameter RI as an edge of an entrainment plateau (phase 
boundary) is approached. As we shall show, the scaling has 
a universal form dictated by the type of bifurcation through 
which the fixed points are lost. 

We first analyze the phase transitions and associated scal- 
ing behaviors that occur in the parameter range when the 
return map F is continuous, as in Fig. |3]^a). Upon varying 
RI the fixed point is here lost through a tangent bifurcation. 
The generic behavior near such a phase boundary is modeled 
by the simple map g(x) = x — x'^ + X, which has a stable fixed 
point at X* — y/\ that is lost through a tangent bifurcation 
as the control parameter A ^ 0. Let 6xn — Xn — x*; then 
5xn+i — Sxn — — 2\/A 8xn — ((5a;„)^, which has large n solu- 
tion &Xn ^ exp(— 2\/A n) * with a coherence time 



3 



^ ~ 1/%/A. This generic scaling holds for each fixed point of 
F'^it) —pTdrv and as the phase boundary is approached from 
within a plateau the coherence time then scales as 



Cr(i?/) 



\RI 



1 



(8) 



consistent with classical exponent = ^ in equilibrium crit- 
ical phenomena [B]. 

The coherence time diverges at the phase boundary (when 
the control parameter RI ~ RI^), and the dynamics can 
again be modeled by the map g{x). In this case A = 0, 
the fixed point x* — and the Sx'^ term above becomes 
relevant. The large n solution is now Xn ~ 1/n and hence 
Xn+i — Xn ^ —1/"-^; in particular the fixed point at the tan- 
gent bifurcation is no longer approached geometrically. Anal- 
ogously, the p/g-entrainment coherence at the phase bound- 
ary then develops according to the power law 



1 



(9) 



consistent with a critical exponent ?/ = 0. 

As we vary the control parameter RI so as to exit the p/q 
entrainment plateau, "bottlenecks" develop near the loca- 
tions of the q lost fixed points. The resulting dynamics can 
again be modeled by the map g{x) which has a bottleneck 
near x — Q for small negative A. As A ^ 0~, the num- 
ber of iterations Nx needed to pass through a fixed interval 
[— c, c] around zero scales like Nx ~ l/-y/|A|. A^a character- 
izes how long in takes to pass through a single bottleneck 
and introduces a time-scale outside the entrainment plateau 
that diverges similar to the coherence time in eqn. (|8|. 

The approach to p/q entrainment is illustrated in Fig. |4] 
which demonstrates exponentially fast coherence of the form 
([t]) inside the 1-1-plateau, the power-law form ([9| at the 
phase boundary and bottleneck behavior just outside the 1- 
1-entrainment phase. The repeated bottleneck behavior rep- 
resents failed entrainment to this particular p/g-entrainment 
phase. If this pattern eventually repeats periodically, the av- 
erage period will converge to different rational multiple of 
Tdrv and hence lie on a different p'/g'-entrainment plateau; 
otherwise Tave/Tdrv is irrational. 

For RI just outside the plateau, the deviation of the aver- 
age period from rTdrv is inversely proportional to the number 
of iterations required to pass through the bottleneck caused 
by lost fixed point of — pTdrv', thus 



\RI-RI^\^. 



(10) 



This average deviation of entrainment plays the role of a 
disorder parameter which, in analogy to equilibrium phase 
transitions, identifies the exponent /3 = ^. This scaling form 
is illustrated in Fig. [5j 

Generically, a smooth map F can acquire or lose periodic 
points only through tangent bifurcations, which dictate the 
universal scaling near the phase boundaries described above. 
However, when the map is discontinuous it can also acquire or 
lose periodic points through a discontinuous bifurcation. As 
we shall see, the associated scaling near the phase boundaries 
for these bifurcations belongs to a new universality class. 
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FIG. 4: Approach to entrainment in ISPs (A^'^) inside (green), 
at (blue) and outside (red) the left edge of the 1/1 plateau, re- 
spectively exhibiting scaling of Q, ([9]l and bottleneck behavior. 

An example of this type of bifurcation is illustrated in 
Fig. [3][b), for r = 2. The behavior at the bifurcation (phase 
boundary) depicted here differs fundamentally from the con- 
tinuous case in Fig. |3ja), in that the slope of the map at the 
fixed point here remains strictly less than one. Consequently, 
iterates of the map still converge to the fixed point following 
the geometric form in ^ , corresponding to a finite coherence 
time at the bifurcation, as opposed to power law scaling. 

To explore the behavior near the plateau edges for dis- 
continuous bifurcations, we introduce a simple map h of the 
form h(x) = ax + X for a; > (with < a < 1) and with 
a discontinuous jump at x = 0. As A sweeps through zero, 
the fixed point of h(x) is lost in a bifurcation similar to that 
seen in Fig. [3|b). Since the map is linear for x > 0, suc- 
cessive iterations satisfy Xn — h^{xo) = x* + a"(a;o — x*) 
while Xn > 0, where x* — A/(l — a) which is a fixed point 
of ft. (x) for A > 0. A bottleneck near x — again develops 
for A small negative, and even though x* < is then not 
a fixed point of h(x), it still controls the passage through 
the bottleneck in terms of the expression for Xn above. The 
number of iterations Nx needed to pass through an interval 
[0, c] is determined by solving — h^^{c) for n (and rounding 
up to the nearest integer). As A ^ 0^, the solution scales 
like Nx ln|A|. 

Consequently, for RI just outside the edge of a plateau 
at which a discontinuous bifurcation occurs, the deviation of 
the average period from rTdrv scales as 



ir„. 



rF 



-1 



In - RI? 



(11) 



So our disorder parameter vanishes logarithmically at phase 
boundaries determined by discontinuous bifurcations, slower 
than any power law. 

Note that the functions F graphed in Fig. [sjjb) are increas- 
ing and have slope +oo on the left at the discontinuities. 
Furthermore, the dependence of F on the parameter RI is 
such that as RI increases, the graph of F moves downward. 
These properties hold in general for all maps F maps under 
consideration as well as the iterates F*. At a bifurcation 
the map F'^(t) — pTdrv still has fixed points but must also 
lie completely on one side of the line y = t. The bifurcation 
is discontinuous if these fixed points occur at the disconti- 
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FIG. 5: Scaling in Tave near the edges of the r = 2 and r — 1 
plateaux. Edge E2 demonstrates the logarithmic scaling in (111 
at a discontinuous bifurcation, while edges E2, and E^ exhibit 
power law scaling in ( 10 1 at tangent bifurcations. 



nuities of the map. Since has slope +00 on the left at 
the jump discontinuities, the fixed points at a discontinuous 
bifurcation must occur on the right side of the jump discon- 
tinuities, and the map F'^(t) ~ pT^rv niust lie below the line 
y = t. These fixed points are lost upon increasing RI and 
consequently discontinuous bifurcations can only occur on 
the right edges of the plateaux, as is the case in the particu- 
lar example illustrated above for r — 2. In other words, the 
bifurcations occurring at the left edges of the entrainment 
plateaux are all tangent bifurcations, even in the parameter 
range where F is discontinuous. 

On the other hand, both types of bifurcations occur at the 
right edges in the parameter range where F is discontinuous, 
although tangent bifurcations are quite rare. For example, 
in Fig. |2] tangent bifurcations occur at the right edges of 
the plateaux for r = |, | and ^; all the other right-edge 
bifurcations we investigated in this parameter range are dis- 
continuous. Furthermore, in this example, the bifurcations 
at the right edges are all discontinuous for RI below some 
threshold. In fact, we can show that such a threshold exists 
if the oscillatory drive E is sufficiently small relative to A?; 
(technically when E < Av/{1 + ((wr)^ -|- l)"^)). The point 
is that under this condition, the map F is concave up every- 



where for RI sufficiently small and since F is increasing, the 
same holds for all its iterates. This rules out the possibility 
of tangent bifurcations at right edges of plateaux; i.e., where 
the map F'^it) — pTdrv is below the line y = t. 

In conclusion, periodically driven IF neurons lose p-q 
entrainment through two markedly different routes, corre- 
sponding to the tangent and discontinuous bifurcations de- 
scribed above. Each has its own characteristic universal scal- 
ing laws which measure the rate of convergence to entrain- 
ment within the p-q plateau as well as the deviation from p-q 
entrainment just outside the plateau. 

In equilibrium statistical mechanics, a phase transition is 
classified as either continuous, which has a diverging coher- 
ence scale (^) and a vanishing order parameter at its critical 
point, or first order, which has no diverging coherence scale 
and a discontinuous jump in the order parameter at its tran- 
sition. The scaling laws at tangent bifurcations are identical 
to that of a continuous phase transition with 'classical' ex- 



ponents j3 



1 

2' 



\ and rj — Q. However, the behavior 



at the discontinuous bifurcations does not match our con- 
ventional understanding of phase transitions with universal 
scaling. The finite coherence time at the discontinuous bifur- 
cations is a feature of first order phase transitions, which have 
discontinuous jumps in their order parameters and no uni- 
versal scaling laws. But as we have seen, Tave varies continu- 
ously and hence the disorder parameter Tave — I'Tdrv vanishes 
at the bifurcation, as is the case for continuous phase transi- 
tions. The logarithmic scaling law for Tave ~ fTdrv vanishes 
more slowly than any power law, and hence exhibits behavior 
which is intermediate between conventional continuous and 
first order phase transitions. 

For conventional continuous phase transitions, universal 
scaling results from the singularity associated with a diverg- 
ing coherence scale. Our discontinuous bifurcation does not 
have such a diverging scale, yet exhibits a new kind of uni- 
versal scaling. Here it is the singularity in the discontinuous 
map that is responsible for universality. 

From a neuro-physiological perspective, our results suggest 
that when a neuron firing with some rate receives an addi- 
tional rhythmic stimulus, a range of p-q entrainment possi- 
bilities exist. Moreover, the convergence to entrainment to 
a p-q phase-locked firing pattern is characterized by a coher- 
ence time which depends sensitively on the distance to the 
entrainment plateau edge. Whole-cell slice recording, where 
an individual cell is stimulated with a constant plus oscilla- 
tory current injection, would be an ideal setting in which to 
explore in vitro the scaling and pattern formation discussed 
in this letter. 
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